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We study the population dynamics of a Bose-Einstein condensate in a double-well potential 
throughout the crossover from Josephson dynamics to hydrodynamics. At barriers higher than 
the chemical potential, we observe slow oscillations well described by a Josephson model. In the 
limit of low barriers, the fundamental frequency agrees with a simple hydrodynamic model, but 
we also observe a second, higher frequency. A full numerical simulation of the Gross-Pitaevskii 
equation giving the frequencies and amplitudes of the observed modes between these two limits is 
compared to the data and is used to understand the origin of the higher mode. Implications for 
trapped matter- wave interferometers are discussed. 
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Quantum mechanical transport is a consequence of 
spatial variations in phase. Supernuids behave like per- 
fect inviscid irrotational fluids, whose velocity is the gra- 
dient of a local phase, so long as the confining potential 
is smooth on the scale of the healing length. Where the 
density is small, as it is near surfaces, quantum kinetic 
terms must be added to the classical hydrodynamic equa- 
tions. Macroscopic quantum coherence phenomena, such 
as Josephson effects, emerge when superfluids are weakly 
linked across such a barrier region. 

Josephson effects have been demonstrated with super- 
conductors pp, liquid helium [2, 3 , and ultracold gases 
in both double- well [U [5] and multiple- well optical trap- 
ping potentials [6]. The canonical description of these 
experiments employs a two-mode model [7H9] , in which a 
sinusoidal current-phase relationship emerges. Hydrody- 
namics has also been studied in both liquids and ultra- 
cold gases [10]. The relative diluteness of gases makes a 
satisfying ab initio description possible [TT] . 

In this Letter, we study the transport of a Bose- 
Einstein condensate (BEC) between two wells separated 
by a tunable barrier and observe the crossover from 
hydrodynamic to Josephson transport. As the barrier 
height Vh is adjusted from below to above the BEC chem- 
ical potential, /i, the density in the link region decreases 
until it classically vanishes when Vb = /i- The healing 
length in the link region, £, increases with and dic- 
tates the nature of transport through this region. Os- 
cillatory dynamics spanning three octaves are observed 
as we smoothly tune £ from 0.3d to 2d, where d is the 
separation between the wells. 

Examination of the dynamics of an elongated BEC in 
a double well is timely. Recent experiments have created 
squeezed and entangled states by adiabatically splitting 
a BEC [T2UT4] . The degree of squeezing inferred in the 
elongated case [I2j [13] seems to exceed what would be 
expected in thermal equilibrium [14] , raising the possibil- 




FIG. 1: (a) Schematic of atom chip double- well trap. Central 
'Z' wire [18] carries static trapping current, Is = 2 A, which, 
with uniform external fields B ex t = (2.2,0.11,0) mT, results 
in an Ioffe-Pritchard style trap with harmonic trapping fre- 
quencies (uj xo ,zo,ujy ) = 2tt x (1300,10) Hz. Side wires are 
1.58 mm from trap center and carry RF currents with am- 
plitude Jrf- This RF current produces a ^-polarized field at 
the trap location with amplitude Brf = 23.6 ± 0.6 \iT (peak 
Rabi frequency Q = 2tt x (82 ± 2 kHz)). A levitation beam 
(pink) is positioned to provide a force cancelling gravity (z- 
direction) while compressing the sample along y. Atoms are 
trapped 190 |^m from the chip surface, (b) A schematic one- 
dimensional cut at t = —0.5 ms through trapping potential 
along x (solid line) in the presence of linear bias (dashed line) 
and (c) balanced potential at t = 0, with Zq ^ 0. 



ity that out-of-equilibrium dynamics may be important. 
With much remaining to be explored in these systems, 
this work represents the first study of the dynamics in 
the crossover regime. 

Our experiment begins as 87 Rb atoms in the \F = 
2,m^ = 2) ground state are trapped on an atom chip 
and evaporatively cooled in a static magnetic potential 
Bs(r), as described elsewhere [15 . To prevent gravita- 
tional sag and to compress the trap in the weak direction 
(with characteristic trap frequency uo y = 2tt x 95 Hz), 
we add an attractive optical potential with a 1064 nm 
beam. We dress the static potential with an oscillat- 
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ing radio- frequency (RF) magnetic field p~6j[T7] radiating 
from two parallel wires on the atom chip (Fig. [TJa)). In 
the rotating- wave approximation (RWA), the adiabatic 
potential created by the combination of the static chip 
trap, the RF dressing, and the optical force is 



U(r) = m' F sgn(g F )hJ5(r) 2 + ^(r 



1 2 2 



(1) 



where m' F = 2 is the effective magnetic quantum num- 
ber, S(r) = cjrf — \hb9fBs( y ) /h\ is the detuning, 
= |/^B#F^RF,±(r)/2ft| is the RF Rabi frequency, 
B RF ,±(r) = |B s (r) x B RF (r)|/|B s (r)| is the amplitude 
of the RF field locally perpendicular to Bs(r), pb is the 
Bohr magneton, gp is the Lande g-factor, % is the reduced 
Planck's constant and m is the atomic mass. By assum- 
ing the individual wells are harmonic near each minimum, 
calculations show that uj z = 2tt x 425 Hz, and uj x varies 
from 2tt x 350 Hz to 2ir x 770 Hz as we tune from low 
to high barriers. For comparison between theory and ex- 



periment, we account for small corrections to Eq. (S26) 
beyond the RWA [19j[20]. 

After turning on the dressing field at a frequency 
urf = 27r x 765 kHz, where the trap is a single well, we 
evaporatively cool to produce a BEC with no discernible 
thermal fraction. In 20 ms, we adiabatically increase ujrf 
to a new value characterized by So = S(r = 0), such that 
the barrier Vh rises and the dressed state potential splits 
along the x-direction into two elongated traps [2Tj . 

Using a second 1064 nm beam weakly focussed off- 
center in x, an approximately linear potential is added 
across the double- well junction to bias the population 
towards one well (Fig. [ljb)) . By applying the bias 
beam before and during the splitting process, we pre- 
pare systems of atoms with a population imbalance Z = 
(N R - N L )/(N R + JV L ), where N R (N L ) is the number of 
atoms in the right (left) well. The range of initial popu- 
lation imbalances Zq = Z(t — 0) we use is 0.05 to 0.10, 
small enough to avoid self-trapping [4]. To initiate the 
dynamics, the power of the bias beam is ramped off in 0.5 
ms (faster than the population dynamics) and the out- 
of-equilibrium system is allowed to evolve for a variable 
time t in the symmetric double- well (Fig. [TJc)). 

To measure the time-dependent population Z(t), we 
freeze dynamics by rapidly increasing both B R y and cjrf 
to separate the wells by 70 |Ltm, where Vh/p ~ 10 4 . We 
release the clouds from the trap and perform standard 
absorption imaging along y after 1.3 ms time-of-flight 
(Fig. J2jb)). Analysis of these images allows us to deter- 
mine A/r and A^l to a precision of ±50 atoms. 

Upon release of the potential bias, we find that the 
population Z(t) oscillates about Z = (Fig. [2|a)) 
[22] . To analyze the dynamics, we use a Fourier trans- 
form (FT) to find the dominant frequency components 
(Fig. |2^c)). We repeat this measurement at many values 
of Vb//i, where p is the Thomas-Fermi chemical poten- 
tial, by varying So- For the purposes of this analysis, we 
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FIG. 2: (a) Population imbalance, Z, vs. time for 5 = 
2tt x (0.1 ± 0.5) kHz, N = 5900 ± 150. The dashed line is 
a decaying two-frequency sinusoidal fit to the data, using two 
fixed frequencies from the FT (lower inset). Each point is the 
average of six repetitions of the experiment; error bars are 
statistical, (b) Averaged absorption image after separation 
and 1.3 ms time of flight, with right and left measurement re- 
gions (dashed boxes) indicated, (c) FT amplitude spectrum 
of data showing two distinct peaks at 268 =b 6 and 151 =b 13 
Hz rising above the noise floor (grey). 



ignore the decay of this signal, the 1/e time constant of 
which is typically two oscillation periods. 

When the barrier is low, Z(t) consistently displays two 
dominant frequency components. For higher barriers, the 
amplitude of the higher-frequency mode decreases un- 
til only a single frequency rises above the noise floor. 
The white points in Fig. [3] give these frequencies as a 
function of the experimental parameter 5o and the calcu- 
lated ratio of barrier height to chemical potential, Vh/ p. 
The ensembles used in Fig. [3] had total atom number 
N = 6600 ± 400 (±1700), where the error bar is statisti- 
cal (systematic). 

In the low- and high-barrier limits, simple models can 
be used to understand the dynamics. For low barriers, 
the hydrodynamic equations of motion can be used to 
estimate the frequency of population oscillation. Assum- 
ing a harmonic population response for some Zq, the re- 
sponse frequency is 
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where p is the density of the condensate at t = 0, S is 
the surface in the y-z plane bisecting the double well, 
and n is the vector normal to this surface. Plotting cjrd 
in Fig. [3] (dotted line), we find good agreement with the 
lower frequency mode at low barriers. Since tunnelling 
cannot contribute to hydrodynamic transport, cjhd ~^ 
as Vb —> p>- The breakdown in hydrodynamics also 
coincides with an increasing healing length, as shown in 
the inset of Fig. [3] 

In the opposite limit, when tunnelling dominates trans- 
port, a Josephson model [8] accurately predicts the fre- 
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FIG. 3: Frequency components of population imbalance vs. 
RF detuning (measured) and barrier height to chemical po- 
tential ratio (calculated). Experimental points (white circles) 
represent the two dominant Fourier components at each de- 
tuning; error bars represent uncertainty contributed by noise 
in the FT from a single time series, but do not include shot-to- 
shot fluctuations. The spectral weight is represented through 
the color map, which has been linearly smoothed between dis- 
crete values of Vb/M an d darker colors indicate greater spectral 
weight. All calculations use N = 8000 and Zq = 0.075, and a 
single-parameter fit of the data to the GPE curves shifts all 
experimental points by £ s hift = 27T x 5.1 kHz [20] to compen- 
sate for a systematic unknown in Bs(0). Statistical vertical 
error bars are shown, while a typical horizontal statistical 
error bar is shown at Vb//i ~ 0.5. Dashed lines represent 
3D GPE frequencies, the solid line the plasma oscillation fre- 
quency predicted by the Josephson model, and the dotted 
line the hydrodynamic approximation, cjhd/27t. White bars 
at Vh//Ji ~ 0.1 indicate the bounds of the GPE simulation 
corresponding to the systematic plus statistical uncertainty 
in atom number. Inset: ratio of healing length, £, to interwell 
distance, d, as a function of Vh/fi. £ is calculated at the center 
of the barrier. 



quency of the highest barrier points, 

"1=^ e (* e+n W)' (3) 

where AE is the energy difference between the symmetric 
and antisymmetric ground states of the double-well po- 
tential, N is total atom number, and /j,\ oc is the chemical 
potential on one side of the well [20]. The agreement is 
surprisingly good even for Vh just above /i, beyond which 
the frequency decreases exponentially. To our knowledge, 
this constitutes the first direct observation of tunneling 
transport of neutral atoms through a magnetic barrier, 
only inferred, for instance, in Refs. p~2j[23]. 

To explain the crossover behavior and the existence 
of the higher- frequency mode, we turn to numerical 
solutions of a time-dependent three-dimensional Gross- 
Pitaevskii equation (GPE) [3 [24], which should describe 
all mean- field dynamics at T = 0. The slope and separa- 



tion of the measured frequencies are well captured by the 
GPE, as shown in Fig. [3J though the decay of population 
imbalance is not reproduced by these simulations. 

The structure and origin of the higher-lying dynamical 
mode can be studied within the simulations. If our trap 
were smoothly deformed to a spherical harmonic poten- 
tial, the two observed modes would connect to odd-parity 
modes [11] : the lower mode connects to the lowest m = 
mode (coming from the t = 1 mode at spherical symme- 
try, where the quantum numbers t and m label the angu- 
lar momentum of the excitation and its projection along 
the axis of symmetry, y, respectively), while the higher 
mode originates from the lowest m = 2 mode (£ = 3 at 
spherical symmetry) [25 . 

With insight from GPE simulations, the observation 
of a second dynamical mode, which was not seen in pre- 
vious experimental work [4, 5 , can be explained. In a 
purely harmonic trap, a linear bias excites only a dipole 
mode [26]. By breaking harmonicity along the splitting 
direction, x, the barrier allows the linear perturbation 
(JL = l,m = 0, where x is the azimuthal axis) to excite 
multiple Bogoliubov modes [27] . Numerical studies show 
that two additional ingredients are required to excite the 
higher mode. First, atom- atom interactions couple the 
x-excitation to the transverse (y, z) motion through the 
nonlinear term in the GPE. Second, the anisotropy of the 
trap in the y-z plane mixes the m = and m = 2 modes 
such that each of the resulting modes drives population 
transfer between wells. 

Figure [4] shows the relative strength R\ = a\j[a\ -fo^) 
of the lower frequency mode as a function of the barrier 
height. The amplitude a\ (0,2) of the lower (higher) fre- 
quency mode is extracted from a decaying two-frequency 
sinusoidal fit. The modes have comparable strength, even 
in the linear perturbation regime, when the barrier is 
below the chemical potential. The small spread in the 
GPE amplitudes shown by the grey band indicates that 
the higher mode is excited independently of the initial 
imbalance, and is not simply due to a high- amplitude 
nonlinearity. 

The trend in R\ reflects the shape of the trap. When 
the barrier is raised from zero, the higher mode is at first 
more easily excited due to an increased anharmonicity 
along x as the trap bottom becomes flatter. By further 
increasing the barrier, the higher- frequency mode disap- 
pears from the population oscillation spectrum due to the 
vanishing excitation of transverse modes. As the wave- 
functions in each individual well are increasingly localized 
to the effectively harmonic minima, the linear bias no 
longer excites intrawell transverse motion. Furthermore, 
in the linear perturbation regime, the interwell Joseph- 
son plasma oscillation, like all Bogoliubov modes, cannot 
itself trigger any other collective mode. 

In conclusion, we have studied the quantum trans- 
port of a BEC in a double-well potential throughout 
the crossover from hydrodynamic to Josephson regimes. 
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FIG. 4: Fraction of low- frequency mode in population dy- 
namics. Dashed line shows the GPE simulation for 8000 
atoms with initial imbalance z(0) — 0.075. The grey shaded 
area represents the variation of the GPE calculations over the 
range of z(0) = 0.05 to 0.10. The vertical error bars are sta- 
tistical; the statistical uncertainty in 6 is 2tt x 0.5 kHz (not 
shown). The GPE calculation gives Ri = 1 when Vb//i — 1.1. 



Apart from fundamental interest, knowing and control- 
ling the nature of superfluid transport is crucial for tech- 
nological applications of weak-link based devices, such 
as double-slit interferometers [121 \21\ I28H30] . The adia- 
batic transformation of a BEC from a single- to a double- 
well trapping potential has been discussed in recent ex- 
perimental works p~2j [14j |23j [3TJ [32] in the context of 
the Josephson model, valid at high barriers [33]. Our 
work demonstrates that for Vb < /i, the lowest mode fre- 
quency will lie below that estimated by the Josephson 
model. Furthermore, the higher-lying mode we observe 
approaches the lowest collective mode as uo y <C oj z [20] 
and may be important to the dynamics of splitting in 
strongly anisotropic double wells [12] [34] . Whether us- 
ing splitting to prepare entangled states [14], or recom- 
bination [32 to perform closed-loop interferometry [31], 
an improved understanding of double- well dynamics pro- 
vides a foundation for controlling mesoscopic superfluids. 
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the GPE. This work has been generously supported by 
ClfAR, CFI, CQIQC, and NSERC. 



[9] 



[1] P. W. Anderson and J. M. Rowell, Phys. Rev. Lett. 10, 
230 (1963). 

[2] S. Pereverzev, A. Loshak, S. Backhaus, J. Davis, and 



R. Packard, Nature (London) 388, 449 (1997). 

[3] S. Backhaus, S. Pereverzev, A. Loshak, J. Davis, and 
R. Packard, Science 278, 1435 (1997). 

[4] M. Albiez et a/., Phys. Rev. Lett. 95, 010402 (2005). 

[5] S. Levy, E. Lahoud, I. Shomroni, and J. Steinhauer, Na- 
ture (London) 449, 579 (2007). 

[6] F. S. Cataliotti et a/., Science 293, 843 (2001); O. Morsch 
and M. Oberthaler, Rev. Mod. Phys. 78, 179 (2006); see 
also B. P. Anderson and M. A. Kasevich, Science 282, 
1686 (1998). 

[7] B. D. Josephson, Physics Letters 1, 251 (1962). 
[8] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, 
Phys. Rev. Lett. 79, 4950 (1997); I. Zapata, F. Sols, and 
A. J. Leggett, Phys. Rev. A 57, R28 (1998); S. Raghavan, 
A. Smerzi, S. Fantoni, and S. R. Shenoy, Phys. Rev. A 59, 
620 (1999); D. Ananikian and T. Bergeman, Phys. Rev. A 
73, 013604 (2006). 

A. Smerzi and A. Trombettoni, Phys. Rev. A 68, 023613 
(2003). 

[10] D. S. Jin, J. R. Ensher, M. R. Matthews, C. E. Wieman, 

and E. A. Cornell, Phys. Rev. Lett. 77, 420 (1996); K. M. 

O'Hara et a/., Science 298, 2179 (2002); T. Bourdel et a/., 

Phys. Rev. Lett. 93, 050401 (2004). 
[11] S. Stringari, Phys. Rev. Lett. 77, 2360 (1996); D. A. W. 

Hutchinson and E. Zaremba, Phys. Rev. A 57, 1280 

(1998). 

[12] G. B. Jo et al, Phys. Rev. Lett. 98, 030407 (2007). 

[13] G.-B. Jo et a/., Phys. Rev. Lett. 99, 240406 (2007). 

[14] J. Esteve, C. Gross, A. Weller, S. Giovanazzi, and M. K. 
Oberthaler, Nature 455, 1216 (2008). 

[15] S. Aubin et a/., Nature Phys. 2, 384 (2006). 

[16] Y. Colombe et al, Europhys. Lett. 67, 593 (2004). 

[17] I. Lesanovsky et a/., Phys. Rev. A 73, 033619 (2006). 

[18] J. Reichel, App. Phys. B 74, 469 (2002). 

[19] S. Hofferberth, B. Fischer, T. Schumm, J. Schmiedmayer, 
and I. Lesanovsky, Phys. Rev. A 76, 013401 (2007). 

[20] See supplementary information. 

[21] T. Schumm et a/., Nature Phys. 1, 57 (2005). 

[22] When the average value of Z differs from zero, we sub- 
tract the average Z from all values of Z{t). In all experi- 
ments, \Z\ < 0.05. 

[23] K. Maussang et al, Phys. Rev. Lett. 105, 080403 (2010). 

[24] L. P. Pitaevskii, Sov. Phys. JETP 13, 451 (1961); 
E. P. Gross, Nuovo Cimento 20, 454 (1961); J. Math. 
Phys. (N.Y.) 4, 195 (1963). 

[25] We checked this numerically by deforming our trap into 
a fully harmonic axially symmetric trap, and following the 
mode frequencies throughout this process. 

[26] W. Kohn, Phys. Rev. 123, 1242 (1961). 

[27] H. Ott et a/., Phys. Rev. Lett. 91, 040402 (2003). 

[28] Y. Shin et a/., Phys. Rev. Lett. 95, 170402 (2005). 

[29] C. Gross, T. Zibold, E. Nicklas, J. Esteve, and M. K. 
Oberthaler, Nature 464, 1165 (2010). 

[30] F. Baumgartner et al, arxiv:1008:1252 (2010). 

[31] Y.-J. Wang et a/., Phys. Rev. Lett. 94, 090405 (2005). 

[32] G.-B. Jo et al, Phys. Rev. Lett. 98, 180401 (2007). 

[33] A. J. Leggett and F. Sols, Phys. Rev. Lett. 81, 1344 
(1998). 

[34] A. Polkovnikov and V. Gritsev, Nat. Phys. 4, 477 (2008). 



5 



Supplementary materials for "Dynamics of a tunable superfluid junction" 



JOSEPHSON MODEL 

We compare our experimental results to the Josephson 
model (JM) and its plasma frequency, uj p (Fig. 3). The 
JM employed here is based on the nonlinear two-mode 
ansatz used in [9], 

tt(r, t) = ^M fl (r; N R (t)) + ^Mr; N L (t)) (SI) 



where i/j RiL = y/N RiL (t) exp(i0 RiL (i)) and </> RiL = (0+± 
(j)-)/y/2 is a real function localized in the left (right) 
well, with ((/)-) being the ground (first antisymmet- 
ric) state of the GPE along the splitting direction. The 
linearized equation gives Z(t) = Z(0) cos(u p t + A0(O)), 
where AO = Or — 9l and the plasma frequency 



^ae (ae + n^ 



(S2) 



where 



dr 



(S3) 

and AE = E_ - E + = 2(/C + N\) with 



E± 
K. 
X 



I 



dr 



2m 



(V0±) 2 + t/(r)4 + - 5 iV4 



dr 



/ 



r a 2 



2m 



(V0 fl )(V0 L ) + ^C/(r)0 L 



(S4) 
,(S5) 
(S6) 



The plasma frequency cu p depends on the derivative 
of the single- well chemical potential /ii OC5 an d therefore 
takes into account the effect of transverse degrees of free- 
dom on the effective nonlinearity determining the inter- 
action energy. This provides an important correction to 
the plasma frequency, typically around 20%. Here AE 
is the energy splitting between the ground state and the 
lowest antisymmetric state along the splitting direction. 
In our experiments, for example, p\ oc /h = 27r x 1.8 kHz 
and AE/h = 2tt x 3.7 Hz at S = 6.9 kHz. 



HYDRODYNAMIC MODEL 

To determine the behaviour of the condensate in a dou- 
ble well in the hydrodynamic regime, we use the conti- 
nuity equation and the equation of motion for the con- 
densate in the hydrodynamic description, ignoring the 



quantum pressure term: 
dp(r,t) 



m 



dt 
<9v s (r,f) 
dt 



V-[v B (r,*)p(r,*)] = 



(S7) 



V [U(r)+gp(r,t) + ^mv s 2 (r,t)] = 



(S8) 



where p(r, t) is the local density and v s is the superfluid 
velocity. We assume harmonic motion of the population 
balance between the wells such that 



^HD 



z 



(S9) 



where cjhd is the hydrodynamic frequency that charac- 
terizes the system. 

The first time derivative of Z = 2N-&/N is 
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TV 



pd 3 r 
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Jv R 

f V-(pv a ) d 3 r = - [ h'(pv a ) dS 
Jv R Js 



(S10) 



where V R is the volume of the right well, S is the area 
of the plane separating the two wells, and n is the unit 
normal vector for this plane. 

The second derivative of Z is then 



Z 



pv s ) dS. 



(Sll) 



To evaluate the frequency, co>hd, we assume that the sys- 
tem begins at rest, such that v s (t = 0) = 0. The time 
derivative of v s is given by the hydrodynamic equation 



of motion, Eq. (|S8j), and 
Z 



t=o 
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mN 



[ pn-V(U(r)+gp) dS. (S12) 
Js 



The geometry of this double well system is such that the 
normal vector ft = x, and the only component of the gra- 
dient which contributes is the x-component. Assuming 
some initial imbalance, Zq, the frequency with which the 
populations oscillate is given by 
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o I Is 



d 



P -K- (U(r)+gp) dydz. 
s ux 

(S13) 

We calculate this initial density profile in the trap, 
tilted by a linear bias Gx, using the Thomas-Fermi ap- 
proach: 



PTF(r) = (fjL-(U(r) + Gx))/g, 



(S14) 



The gradient term in the integrand of Eq. ( |S13| ) is then 
simply — G. 



6 



The characteristic frequency is thus 
2G 



^HD 



mNZo 



p TF dydz, 



(S15) 



which indicates that the frequency can be found by sim- 
ply evaluating the density at the surface between the 
two wells and integrating over the region by which the 
two halves are connected. From this expression, we see 
that the uj^ d decreases as the area connecting the wells 
decreases, and falls to zero when the barrier surpasses 
the chemical potential and the Thomas-Fermi density is 
strictly zero on the plane S. 



The equation (S13), upon substituting ptf with the 



Gross-Pitaevskii ground state density, is also valid when 



we include a quantum pressure term in Eq. (S8). How- 
ever, even with the quantum pressure, this model is not in 
exact agreement with the Gross-Pitaevskii equation, due 
to the non-harmonic component in the oscillation. At 



high barriers, though anharmonicity is small, Eq. (S13) 
is less accurate than the JM. 

We have checked that the frequencies predicted by 



Eq. (S15) are consistent with the dynamical simulations 



of Eqs. (S7), (S8), done using a test particle method [S2]. 



GROSS-PITAEVSKII EQUATION 

We solve numerically the time-dependent equation 

h 2 



2m 



V 2 *(r)+ (Z7(r)+#|*(r)| 2 )*(r) 

(S16) 

where \£(r) is the complex condensate order parameter, 
U (r) is the double- well external trapping potential, and 
g = 4:7rh 2 a s /m with a s the 87 Rb s-wave scattering length. 
While all calculations were done using Eq. (S26), an in- 
tuitive understanding of the potential emerges from the 
separable approximate form: 



Us. 



1 2 2 

^muj y y 



2^2 



a 2 x 



(S17) 



where a 2 < 0. The trap exhibits an axial anisotropy, with 
elongation along the ^/-direction such that uo z ~ Auo y . 



DECAY OF POPULATION IMBALANCE 

In comparing the time series measured in the exper- 
iment with those found from GPE calculations, we see 
similar multiple- frequency behaviour. One striking dif- 
ference is the presence of "decay" in the experimental 
data - the fall off of the amplitude of the populations os- 
cillations with time. The characteristic time scale of the 
decay, r, is approximately equal to two oscillation periods 
over all values of Vb/ja. We model this as an exponen- 
tially decaying envelope in our analysis, and include it in 
our fitting equation (Eq. ( |S27| )). 



N 

CD 
O 

a 

x> 
B 

a 

o 

I — I 

Oh 
O 



0.1 
0.05 


0.05 



-0.1 



1 1 


1 1 


f\li 


1 A- 













5 10 15 20 
hold time (ms) 

FIG. SI: Comparison of experimental and GPE time series for 
8 = -0.1 kHz (GPE)) and S = -0.1 ± 0.5 kHz (experiment). 
Experimental points are shown as black dots, and the fit to 
experimental data is shown as a black dashed line. The GPE 
results are shown as a solid blue line. 



In the GPE results, no such decay is observed. Figure 
[ST] shows a comparison between one experimental run 
and a GPE calculation for very similar parameters (6 = 
—0.1 kHz). Indeed, GPE calculations to 64 ms show no 
sign of damping. Besides the possibility of the damping 
arising from technical sources, it may be due to thermal 
or other stochastic effects not included in the T = mean 
field calculation. 



ROLE OF TRAP ANISOTROPY 



We studied the role of the trap anisotropy (i.e., uj z ^ 
uj y ) by observing the transformation of the m = and 
m = 2 modes as the trap is deformed from axially sym- 
metric to strongly axially anisotropic, in presence of a 
purely anharmonic potential along x, using the simpli- 
fied potential, Eq. (S17). The dipole perturbation excites 



both modes as soon as the axial symmetry is broken, 
and the spectrum shows a second frequency growing in 
strength as the axial anisotropy is increased. The values 
of the mode frequencies as a function of uo z /uj y are shown 
in Fig. [S2| Close to axial symmetry, the lower frequency 
depends only slightly on the transverse confinement, indi- 
cating that the m = mode is a dominant component of 
the Bogoliubov excitation. Motion is primarily along the 
splitting direction without oscillations in the transverse 
directions. 

Sufficiently far from axial symmetry, both frequencies 
start to decrease with increasing anisotropy and show a 
similar behavior. In particular, the experimental trap- 
ping conditions correspond to the point uo z « 4uj y , as 



indicated in Fig. [S2j where the two modes begin to show 
a similar dependence on transverse confinement. This 
strongly suggests that for such high axial anisotropy, each 
Bogoliubov mode is mainly a combination of the two orig- 
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FIG. S2: Mode frequencies for m = (open) and m = 2 
(closed) as a function of trap anisotropy. Grey arrow indicates 
the anisotropy used in this experi ment . These calculations use 
the approximate potential Eq. (S17). These are calculated 
with fixed uj z and decreasing uo y 



inal m — and m = 2 modes at axial symmetry. 



FULL DESCRIPTION OF RWA POTENTIAL 

The double- well potential is created through a coupling 
between static and rf magnetic fields. In the dressed state 
picture, these combine to form the effective potential [T7] 



#FMB#rf,±(r) 



(S18) 

where m' F is the adiabatic magnetic quantum number, gp 
is the Lande factor, \±b is the Bohr magneton, B$ (r) is 
the static magnetic field, described by an Ioffe-Pritchard 
potential, and 7? r f,±(r) = |B s (r) x B rf (r)|/|B s (r)| is the 
component of the oscillating magnetic field locally per- 
pendicular to the static field at each point, r. 

The static magnetic trap arises a result of the com- 
bination of current flowing through the 'Z'-wire on the 
chip, an external bias field, and an external Ioffe field. In 
combination, these create an Ioffe-Pritchard style trap, 
a static magnetic field B$ = B X Z + B y y + B z z, whose 
components are described by 



B x {x,z)=B'x-^xy, (S19) 

B y (x,y,z) = B S (0) + ^(y 2 - + ^)) , (S20) 
and (S21) 

(S22) 



B" 

Bz(y,z) = -B'z- -yyz. 



In the limit of a small cloud, the static potential is well- 
approximated by a harmonic trap, characterized by ra- 
dial and axial trapping frequencies uj XjZ and uj y . In terms 



of these measurable values, the static trap-bottom term, 
the gradient term, and the curvature term are given by 



Bs (0) 



B' 



B h 



2/zcjtb 
rap^FMB 



mB s (0) 
ra F # F /i B 



2 



ra F # F /iB 



(S23) 
(S24) 
(S25) 



respectively, where we define cjtb = Mb#f^s(0)/^ as the 
"trap bottom" frequency. 

Typical values for the parameters in Eqs. ( |S25[ ) and 
( |S26| ) are: (u x ^u y ) = 2tt x (1310,10) Hz, u TB = 2tt x 
787 kHz, £ rf ,± = 240 mG, u yi0 = 2ir x 95 Hz, and in the 
\F = 2, nip = 2) state of 87 Rb we use, m f F g^ = 1. 



CORRECTIONS TO THE ROTATING- WAVE 
APPROXIMATION 



To calculate our trapping potential, Eq. (|S26| as- 
sumes the rotating- wave approximation (RWA), but as 
discussed in [19], the RWA fails for large Rabi frequen- 
cies. We study the effect of the beyond- RWA effects for 
our trap and find that we can account for the difference 
between the approximate and full potentials by simply 
shifting the detuning by a fixed amount. 

We calculate the full potential in a two-dimensional 
plane at y = for our trap at a particular detuning, So, 
as described in [19]. This 2D contour is fit to 



U RWA = m^gn^M^r) - ^shift) 2 + (S26) 



at y = 0, which is just Eq. (S26) without the compression 



term, where 5 s hift is the only fit parameter and describes 
a shift of the detuning. We compare the shape of the full 
potential to the RWA potential with the shift and find 
that they are very similar. The shift is calculated for 
all detunings used in this work and is roughly uniform 
for the range explored (Fig. S3). We apply a shift #o — » 
So + 2tt x 1.9 kHz to each of the detunings used with the 
RWA in this work. 

The shift we find is of the opposite sign to that found in 
Ref. [19 . Compared to the potential used in that work, 
our Rabi frequency is much smaller and our detuning 
much closer to zero, and we have confirmed that the shift 
changes sign for larger Q and large negative detunings. 



SYSTEMATIC TRAP-BOTTOM SHIFTS 

As noted in the main text, we shift the data by a fixed 
detuning for display in Fig. 3. Despite the discrepancy 
between the calculated and measured values of the de- 
tuning, the data matches the GPE in terms of the shape 
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FIG. S3: Comparisons between RWA and full potential cal- 
culations, (a) Potential energy curve through y = z = with 
So/2tt = 0, calculated using full potential (blue dots), RWA 
approximation (red solid line) , and RWA approximation with 
fitted shift of £ s hift/27r = —1.85 kHz (black dashed line). In- 
set: difference between full potential and RWA potential (in 
Hz) with shift over entire 2D plane at y — used for fit. Color 
bar indicates in Hz the difference of potentials, (b) Fitted de- 
tuning shift as a function of detuning, i.e., the number one 
should subtract from the detuning in the RWA expression to 
get the best estimate of the potential. 



of the curves and the slope of the population oscillation 
frequency as a function of Vb/fi. For this reason, we 
are satisfied that the shift we are applying is acting to 
account for an unknown systematic uncertainty in the 
determination of the trap bottom value B$ (0). 

After taking into account all known systematic shifts, 
which include the beyond-RWA effect described above 
and calibrations in measuring the trap bottom B$ (0), 
we fit the experimental data to the GPE simulation data 
using a single-parameter least-squares fit, where the fit- 
ting parameter acts to slide the data along the detuning 
axis. We find that a shift of 2ir x (5.1 ±0.1) kHz accounts 
for the difference between the experiment and the GPE. 
This shift is in the opposite direction to the beyond-RWA 
corrections. Possible sources of this discrepancy include 
systematic errors in determining the static trap bottom 
Bs (0), or imperfections in the polarization of Brf clue 
to the proximity of the fields to the chip and its copper 
support. 



ATOM NUMBER 

To calibrate the atom number, we use standard absorp- 
tion imaging to measure the thermal fraction of clouds 
above and below the condensate temperature. We deter- 
mine the total atom number by measuring the total ab- 
sorption of the cloud, and the thermal number by fitting 
the wings to a Bose-Einstein distribution and integrating 
under the entire curve to extract thermal atom number. 
The temperature of each condensate is determined by 
fitting the wings to a Gaussian. 

To find Tc, the condensation temperature, we plot 



the condensate fraction as a function of temperature. 
We determine the temperature at which the conden- 
sate fraction is first non-zero, and find the number of 
atoms to which this corresponds. Using the relation- 
ship between condensation temperature and atom num- 
ber, including finite size and interaction effects [SI], we 
can determine the condensation temperature to ±9% 
(Tc = 640 ± 40 nK). Propagating this error through to 
atom number, we arrive at a calibration factor iVa CtU ai = 
^measured x (1.3 ±0.3), which accounts for the systematic 
uncertainty in our atom number, N = 6600 ± 1700. 

The number, N = 8000 was chosen for the calculations 
because this is the number within the systematic uncer- 
tainty for which the best agreement is found for mode 
amplitudes (Fig. 4). The same N = 8000 is used in the 
Josephson model and hydrodynamic approximations. 



DATA ANALYSIS 

To analyze the time series data, as in Fig. 2(a), we use 
a Fourier transform (FT). To prepare the data, we elim- 
inate the offset from Z = components by subtracting 
from each point the mean, where the mean might be non- 
zero due to a small equilibrium imbalance in the system. 
We smooth the transformed data by padding the time 
series with zeros to a total of 1024 points. 

When extracting the peak locations from the FT, we 
eliminate the points below the frequency given by 1/ttot? 
where £ t ot is the longest hold time. We then find the two 
peaks with the maximum height and use these are our 
data points in Fig. 3. We plot the amplitudes-squared 
of the FTs in the color-map in Fig. 3 behind the data. 
The data are linearly interpolated numerically between 
the values of V^/fi (So) at which the data were measured. 

The uncertainty in the frequency measurement is found 
by simulating data with the same level of noise as the 
original time series. The quantity of noise is determined 
by fitting the time series to a 2-frequency decaying expo- 
nential function 

Z(t) = e~ t/r [ai sin(27r^i(t - t 01 )) + a 2 sm(27rv 2 (t - £02))] 

(S27) 

where r is a time constant for decay, 0,1(2) is the ampli- 
tude of the first (second) frequency component, and V\{2) 
is the first (second) frequency component, and £01(02) 
is the constant accounting for the phase shift of the 
first (second) component. The standard deviation of the 
residuals from this fit gives the noise level. We simulate 
data 100 times with the same parameters as those given 
by the fit, with the same total time and density of points, 
but with different randomized instances of Gaussian noise 
whose standard deviation is the same as that measured. 
Taking the frequency measurements from each of these 
trials, we determine the smallest range inside of which 
68% of the measurements lie. This confidence interval is 
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used as the uncertainty in the frequency measurement. 

The noise floor in the FT is established in a similar 
fashion. Using the result for the noise level from the 
time series, we simulate pure Gaussian noise and take 
the FT of this. The noise floor we show is the mean plus 
one standard deviation of the maximum peak amplitudes 
found in 100 such simulations. 

The amplitudes used in determining the ratio Ri, 
shown in Fig. 4, are given by the values determined by 



the fit (Eq. (S27)). The uncertainties in these values are 



determined in a similar way to those in the frequencies; 
we use the noise level in the residuals of the fit, simulate 
and fit 100 sets of data with similar parameters, and use 
the 68% confidence interval of these results to represent 
our uncertainty. 

One significant difference between the calculated and 



measured quantities is that the calculated amplitudes dis- 
play no decay. The measured values, which come from 
the fits to the Eq. (S27), rely upon the fitting routine 
to extrapolate backwards in time to deterimine the t = 
amplitudes. The uncertainty associated with this process 
results in the scatter in the measurements, and may be a 
cause of some of the discrepancy between the calculated 
and measured values. 
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